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Abstract 

We analyze the thermalization properties and the vahdity of the Eigenstate Thermahzation 
Hypothesis in a generic class of quantum Hamiltonians where the quench parameter explicitly 
breaks a Z2 symmetry. Natural realizations of such systems are given by random matrices 
expressed in a block form where the terms responsible for the quench dynamics are the off- 
diagonal blocks. Our analysis examines both dense and sparse random matrix realizations of 
the Hamiltonians and the observables. Sparse random matrices may be associated with local 
quantum Hamiltonians and they show a different spread of the observables on the energy 
eigenstates with respect to the dense ones. In particular, the numerical data seems to support 
the existence of rare states, i.e. states where the observables take expectation values which 
are different compared to the typical ones sampled by the micro-canonical distribution. In the 
case of sparse random matrices we also extract the finite size behavior of two different time 
scales associated with the thermalization process. 
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1 Introduction 



Largely triggered by recent experiments on cold atoms [H El O HI |5] , there has been in the past 
few years intense theoretical activity aimed at understanding the non-equilibrium dynamics 
in closed interacting quantum systems following a change in one of the system parameters. 
In many interesting cases, the system parameters can be changed rapidly with respect to 
all other time scales of the system that it is meaningful to consider the limit known as a 
quantum quench: namely, the system is prepared in an energy eigenstate l-i/^o) of an initial 
pre-quench Hamiltonian, Hpre, and then is allowed to evolve according to a new post-quench 
Hamiltonian, Hpost, which differs from -ffp^e by some variation of a parameter. Given that 
such a time evolution is purely unitary, one may wonder whether the system reaches, for large 
times, a new steady state and, moreover, if measurements done on this state are able to be 
related to the thermal density matrix of the system. 

Recent progress in understanding thermalization of an extended quantum system following 
a quench has involved both analytical and numerical studies. From the analytical point of 
view, one of the first results can be traced back to von Neumann [6] (see also 0), who pointed 
out the subtleties involved in defining a quantum generalization of the notion of ergodicity. 
Imagine, for example, that with an appropriate choice of the pre-quench Hamiltonian, Hpre, 
we have prepared the system in a linear superposition of M energy eigenstates, \Ea), of the 
post-quench Hamiltonian Hpost, all in a shell of energies, \Ea — E\ < A, centered at E: 

\Ea-E\<A 

The time average of the density matrix based on this state, given by 

PdiagiE) = W(tWW\ = Yl \Ca\^\Ea){Ea\ , (1.2) 

\Ea-E\<A 

defines the so-called diagonal ensemble which is, in general, different from the micro-canonical 
density matrix defined by 

Pn.c(E) = ^ Y ■ (1-3) 

\Ea-E\<A 

So, unless it happens that |cap = l/A/", quantum ergodicity is strictly speaking almost never 
realized. Things may be however different if one looks at the long time evolution (and time 
average) of the expectation value of macroscopic observables, O. For these quantities, defining 
{0)diag = Tr(C'/9diag(E)) and {0)mc = Tr(C'/9mc(E)), it may be true that the identity 

WWOm)) = {0)d^ag = {0)^mc , (1-4) 

indeed holds. One possibility is that the expectation values Oaa = {Ea\0\Ea) of the macro- 
scopic observables do not fluctuate between the Hamiltonian eigenstates which are close in 
energy. In this case, in fact, the identity (|1.4 p holds for all those initial states which are 
sufficiently narrow in energy. This is, in a nutshell, the scenario known in the literature as the 
Eigenstate Thermalization Hypothesis (ETH)) which was put forward by Deutsch and Sred- 
nicki |8l [9] , based on previous work by Berry |10| , and which has been recently advocated by 
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Rigol et al. as the mechanism behind the thermahzation processes in quantum extended 
systems. 

Recently this hypothesis has been put under intense scrutiny by different groups. The 
main emphasis heretofore has been given to the numerical analysis of specific model^, such 
as hard-core bosons |1H |25] . the Bose-Hubbard model [26], strongly correlated interacting 
fermions |27) . the Hubbard model |28) . one dimensional spin chains [T^, etc. In this paper, 
instead of analyzing a particular system, we take a different approach. Namely, our strategy 
herein is to study the thermahzation properties of a generic class of Hamiltonians and a generic 
set of observables. The natural language to approach the problem from this point of view is 
obviously provided by random matrices jj^, which in the following will parametrize both 
the Hamiltonians and the observable^. In particular we have chosen to study the quantum 
quenches and the relative thermahzation in a class of Hamiltonians given by 

H{h) = Ho + hV , (1.5) 

where the quench parameter h is meant to explicitly break a Z2 symmetry of the 'unperturbed' 
Hamiltonian Hq. Such Hamiltonians, which are arguably among the simplest examples of 
quantum systems, may model spin chains in the presence of an external magnetic field but, 
as we shall see later, they may also encode the familiar quantum Ising chain in a transverse 
magnetic field. Given the relative simplicity of this class of Hamiltonians, studying their 
quench dynamics may be a useful path to extract interesting information on generic properties 
of non-equilibrium systems, thus disregarding, in doing so, all additional complications coming 
from a richer structure of states of a specific model. 

However, even upon adopting the abstract language of random matrices, an important 
issue governing thermahzation properties and of which to be mindful is the locality of the 
interaction. Indeed, the mechanism behind thermahzation may be different in the case in 
which the Hamiltonian is local and when it is not. As argued below, local Hamiltonians 
and local observables are associated to sparse matrices, i.e. matrices with a small proportion 
of non-zero entries, while non-local Hamiltonians and non-local observables are described by 
dense matrices. The two kinds of matrices have two different densities of states and it is 
essentially for this reason that one observes a different behavior of the eigenstates under a 
quench of the parameter h. 

Important features of quantum quench processes in local systems were discussed in a paper 
by Biroli et al. |32) . in particular the role played by rare fluctuations in the thermahzation 
of local observables. These authors considered the existence of certain rare eigenstates - rare 
compared to the typical ones sampled by the micro-canonical distribution - but which may 
be responsible, if properly weighted, for the absence of thermahzation observed in certain 
systems. As discussed in more detail later, the presence of such states can be detected by 
studying the spread of the expectation values of the observables on the energy eigenstates, 

''Analytic results for quantum quenches have been obtained only for a restricted class of exactly solvable 
lattice models, such as the XY chain, the Ising model or the XXZ quantum spin chain |12l 1131 1141 1151 116) . 
Analytic results have been also obtained for systems nearby the critical point |18) or for continuous exactly 
solvable systems, especially in the regime of conformal symmetry |19l 1201 121) . However it has been argued 
that the relaxation phenomena of these models, ruled by an infinite number of conserved quantities, may be 
different from the thermahzation of a generic model and may require the introduction of a generalized Gibbs 
ensemble, as proposed in |22| (see also |23) for a derivation in integrable field theories). In this paper, however, 
we will not deal with such systems, but rather address these issues in a separate publication. 

'^For simplicity we consider hereafter real symmetric matrices. 
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in particular by the finite size dependence of tlie distribution of expectation values. The 
numerical analysis that we have performed seems indeed to indicate the existence of these 
rare states in the case of sparse random matrices, while they are absent in the case of dense 
random matrices. However, in our numerics, thermalization is observed nonetheless in sparse 
matrix Hamiltonians, simply because our averaging procedure on the different sampling of 
observables and Hamiltonians does not place a natural exponentially large weight upon the 
rare states, thus enabling them to break thermalization. 

It should be underlined that the existence of rare states in the thermodynamic limit has 
been debated in the literature and in particular in a series of papers by Santos and Rigol 
|33| IMl |35] . They have argued that in a portion of the phase diagram of an extended t-J 
model with next-nearest-neighbour interactions, rare states are absent. We will come back to 
this conclusion in our presentation of results. 

The paper is organized as follows: in Section |2] we discuss the notion of locality in re- 
gards to a Hamiltonian for a quantum system, and we argue that the corresponding matrix 
representation is given by a sparse matrix. We follow this in Appendix A with an extended 
discussion of whether a matrix sparse in a position basis is also sparse in the corresponding 
momentum basis. In Section |3] we address the issue of the ETH and, following Ref. |32], we 
briefly remind the reader of certain salient issues regarding the application of the ETH to 
our quenches. Section U] is devoted to the discussion of our quench protocol implemented on 
Z2 breaking quantum Hamiltonians. Section |5] deals with this protocol as applied to dense 
random matrices, while Section [6] treats sparse random matrices. In Section [7] we discuss the 
relevant time scales for thermalization in sparse random matrices, while our conclusions are 
presented in Section |8l 

2 Locality 

In this section we discuss the nature of Hamiltonian matrices associated with local models. 
In order to consider finite-size matrices, we will always have in mind that the model contains 
(effectively) proper physical cut-offs (for example, volume or lattice spacings) such that the 
Hilbert space involved is actually finite. For local models here we mean models local in the 
space variable, x, (if continuous) or in the lattice site, i, if discrete. The main idea of this 
section is the following: in a generic basis of the Hilbert space, the matrix representation of a 
local Hamiltonian corresponds to a dense matrix, i.e. a matrix which has all entries different 
from zero (an explicit example will be given below). However, if the theory is local, there 
will exist a basis (in the following called the local basis), in which the Hamiltonian will be 
represented by a sparse matrix, i.e. a matrix where the great majority of its entries are zero. 
In the following, for simplicity, we focus our attention on one-dimensional systems. 

2.1 Local basis 

To clarify the concept of the local basis, it is worth starting our discussion by a paradigmatic 
example: the Id quantum Ising model in a longitudinal field (generically a non-integrable 
model). In this case the quantum Hamiltonian for L sites is given in terms of Pauli matrices 
and takes the form: 

L L 
H^Y. «+i + ^< + < = E (2.1) 

i=l i=l 
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The last equality makes evident the local nature of this model: the Hamiltonian has been 
written as a sum of operators involving only two lattice sites. The matrix representation can 
be obtained once we fix a basis. A typical choice is given by the set of common eigenstates of 
the af operators. They can be written as, 

I tt •••!), in •••;),••• \mim2...mL), (2.2) 

where each rrii £ {t)i} corresponds to the two possible eigenstates of af: 

t) = I t), i) = -\i) =^ (^ilm) = mi\mi). 

Therefore the Hilbert space is made of = 2^ elements. We call this the real-space basis and 
it can be characterized by the fact that its elements are tensor products of the states of each 
site. Now let's consider the matrix element of the Hamiltonian density Hi: 

H.'"" = {mi... mN\Hi\m[ . . . m'^) = {{rmmi+i + mi)5m„m[ + h6^,-m'Jl ^nik,m'^^ (2-3) 

where s, s' are shorthand for the full set of indices, mi, . . . , m^r. From this expression it is 
easy to deduce that on each row of the matrix there are L + 1 non-zero entries and therefore 
the total number of non-zero elements of the N x N matrix H is M = {L + 1)N. Since the 
total number of matrix elements is A^, the density of non-zero elements is given by 

J\f flnN \ 1 InA^ , , 

+ !-«—-, (2.4) 



A2 \ln2 J N N 
while the density of the zeros of the matrix H is 

klnN , , 

po = l-p^l-^^. (2.5) 

The constant k is related to specific properties of the model, but the scaling in Eq. (|2.4 p 
is general. Therefore for large values of A, the Hamiltonian matrix Hj^ is a sparse matrix, 
i.e. a matrix with a large number of zeros and few non-zero entries. This statement holds in 
general for any local quantum Hamiltonian, as for instance the one coming from a quantum 
field theory: the fraction of non-zero elements of the Hamiltonian is exponentially suppressed 
in the thermodynamic limit L — t- oo. Take for instance a (1 -|- 1) bosonic field theory, whose 
Hamiltonian can be written as 



^U\x, t) + ^id,^{x, t)f + V{^{x, t)) 



dx , (2.6) 



where 11 (x, t) = ^ is the canonical conjugate of the operator ip{x,t) entering the equal-time 
commutation relation, 

[ip{x,t),U{y,t)] = i6{x-y) . (2.7) 

Once the space coordinate x and the field values have been discretized, we can choose the basis 
of eigenstates of the operators ifr^. and the Hamiltonian (|2.6 P can be expressed as a matrix. 
All the terms except the momentum become diagonal in this basis. The momentum instead 
plays the same role of a "spin flip" operator as does in the Ising case: one can easily see 
that this term provides off-diagonal coefficients in the matrix, and in each row they will be as 
many non-zero entries as there are lattice sites. Therefore a scaling analogous to Eq. (|2.5 P is 
found. 
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2.2 The basis of momenta 

In the last section we pointed out the kind of matrices one should expect when a local model 
is formulated in a real-space basis. However one can ask the question what happens to the 
sparseness of the matrix if the basis is changed. Among all other possible bases, it is the 
momentum basis that is often used to express the Hamiltonian in matrix form. For this 
reason we analyze here what would be the differences in this case. The details can be found 
in Appendix [Aj where one can see that the sparseness of the Hamiltonian depends on the 
particularities of the momentum basis that one considers. But in summary: 

• If one considers the basis of momenta obtained as a Fourier transform of the rigid 
translation of the real-space basis, the Hamiltonian will appear again as a sparse matrix, 
perhaps even with a larger density of non-zero entries. This is true simply because the 
change of basis we are considering is sparse, i.e. for a state at L = 4: 

\<P) = \ ttu) -^\<Pk) = \ mi) + e"^^/2| im) + e^"'''/^! an) + e"''^i'^\ mt) 

and the number of terms in the last expression is always less than L ^ 2^. 

• If one considers the basis of momenta obtained taking the Fourier transform of the free 
single particle excitations: 

i=o 

the change of basis in each block of defined total momentum will not be sparse. In the 
general case, the Hamiltonian matrix will be characterized by dense blocks of fixed total 
momentum. 

These considerations allow us to stress one point: there does exist a basis, regardless of inter- 
actions, in which momentum is a good quantum number and which maintains the sparseness 
of the Hamiltonian matrix. 

After these general remarks, we now focus on the real-space basis keeping the scaling Eq. 
(12.5 P describing the sparseness of the matrix. 

3 Quantum quenches, thermalization, and the ETH 

Let us consider an initial state l-f/^o) which is an eigenstate of an initial Hamiltonian, H{h^), 
governed by the parameter h^. At t = we abruptly (i.e. non-adiabatically) change the 
value of the parameter to /i^. The evolution of the initial state will be then governed by the 
dynamics given by H{K^). Our interest is in the long time behaviour of expectation values of 
some one-point observable, {tpQ[t)\0\'il)Q{t)) . An observable has a thermal behavior if its long 
time expectation values coincides with the micro-canonical prediction, i.e. 

(Vo(t)IOlV'o(t)) ^Tr(Op^e) = (0)mc . (3.1) 

However as the systems we will consider are finite sized, one has to allow for quantum revivals. 
A proper way to understand the above limit is then to require that eq. (|3.1 P holds in the long 
time limit at almost all times. Mathematically, this means that the mean square difference 
between the LHS and RHS of eq. (13.1 p averaged over long times is vanishingly small for large 
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systems - a restatement of von Neumann's quantum ergodic theorem |6l [7]- We will come 
back to the time scales and their subtleties involved in the approach to equilibrium in Section 
[T] For now we say that an observable thermalizes if its long-time average coincides with the 
micro-canonical average, namely 



where Ca = (tpolEa) are the overlap of the initial state on the eigenstate \Ea) of H{h^), 
and Oaa = {Ea\0\Ea) are the expectation values of the observable, O, on the post-quench 
eigenstates. Eq. (|3.2 p defines the diagonal ensemble prediction, with the corresponding 
density matrix defined as 



supposing the eigenstates of H(h^) are non-degenerate. 

A possible mechanism for the thermal behavior of physical observables is based on the so 
called Eigenstate Thermalization Hypothesis (ETH) jUlH]. It states that the expectation value 
of a physical observable, Oaa = {Ea\0\Ea) , on an eigenstate, \Ea), of the Hamiltonian is a 
smooth function of its energy, E^ , with its value essentially constant on each micro-canonical 
energy shell. In such a scenario, thermalization in the asymptotic limit follows for every initial 
condition sufficiently narrow in energy. ETH implies that thermalization can occur in a closed 
quantum system, different from the classical case where thermalization occurs through the 
interactions with a bath. As pointed out by Biroli et al. |32) there are two possible inter- 
pretation of ETH: a weak one, which can be shown to be verified even for integrable models, 
which states that the fraction of non-thermal states vanishes in the thermodynamic limit, and 
a strong one which states that non-thermal states completely disappear in the thermodynamic 
limit. In the weak version of the ETH, not every initial condition will thermalize. 

We briefly remind the reader of the origin of these two interpretations as it will be salient 
later. Firstly, for thermalization to occur one needs a distribution of the overlaps peaked 
around the energy E = (^ol-f^l^o)- As shown in Ref. [TT], the energy density e has vanishing 
fluctuations in the thermodynamic limit 



where L is the system size and a is the dimension of the space over which the coupling 
h is adjusted in the quench. As our random quench is meant to be a proxy for a generic 
one dimensional quench, a = 1. However we will however see that a is effectively larger 
when we consider quenches in dense matrices, and consequently Ae does not vanish in the 
thermodynamic limit. Correspondingly this property means that the distribution of intensive 
eigenenergies (eigenenergies scaled by 1/L) with weights |cap is peaked for large system sizes. 
If the ETH is true, an immediate consequence of property Eq. (|3.4 p would be that averages 
in the diagonal ensemble coincide with averages in the micro-canonical ensemble. However, 
for a finite system, there will always be finite fluctuations of Oaa- To characterize the ETH 
mechanism we need then to have some control on the evolution of the distribution of Oaa in 




(3.2) 




(3.3) 



a 




(3.4) 
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approaching the thermodynamic hmit. As shown in |32] the width of the distribution Oaa of 

an intensive local observabl^ vanishes in the thermodynamic limit 



where e is the intensive energy defining a micro-canonical shell including \Ea) such that Ea/L G 
[e — e, e + e] and A'^e is the number of states in the microcanonical shell. Eq. (13.5 P implies that 
the fraction of states characterized by a value of Oaa different from the micro-canonical average 
vanishes in the thermodynamic limit. However states with different values of Oaa may exist. 
These states live in the tails of the shrinking Oaa distribution and are expected to be small in 
number. This is why they are called "rare". These states, however, under proper conditions, 
can be relevant to the issue of thermalization. Indeed, if in the |cap distribution they are 
weighted heavily, the diagonal ensemble average will be different from the micro-canonical 
and the system keeps a memory of the initial state. As emphasized in |32) . it is clear that the 
weak interpretation of ETH does not imply thermalization in the thermodynamic limit for 
every initial condition, while, with the proviso that Eq. 3.4 holds, the strong interpretation 
does. 

4 Z2 symmetry breaking quench protocol 

The class of Hamiltonians we choose to study can be thought as akin to the quantum Ising 
model in the presence of an additional longitudinal field. The quench protocol involves, for 
the sake of specificity, taking h to —h. This quench reflects that the Ising Hamiltonian is not 
invariant under the Z2 operator, V = e'^'^i^/'^+Sz)- 



In order to mimic this in the context of random matrices, we suppose we have divided the 
canonical basis of Ising states into two groupings, even and odd under V, and to then have 
sorted them by ordering all even states before any odd states. In Ising, the transverse field 
term couples states with different parity, such that the dependence of the Hamiltonian on the 
external field is seen in the off-diagonal blocks, i.e. 



It is this form then that we take for our random matrices. 

Observables of the systems associated to the Hamiltonian (|4.1 p can be split into even and 
odd Z2 classes. This classification is again motivated by the case of the Ising-spin chain in a 
transverse magnetic field, where the natural observables cj^ and ax are respectively odd and 
even w.r.t. to the action of V. The even observables have non-zero elements in the diagonal 
blocks alone, while the odd observables are non-zero only in the off-diagonal blocks: 




(3.5) 



VH{h)V^ = H{-h). 




(4.1) 




(4.2) 



'^O is an intensive local observables if it can be written as "Yla where Oa are finite ranged observables 
and the sum is over a local spatial region. 
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where the volume factor L has been added to make these quantities intensive. Since we expect 
the Hilbert space to be exponentially large in the volume of the system we fix the system size 
corresponding to an x random matrix via 

L = lniV . 

After defining the Hamiltonian as above, we will analyze the quench dynamics under the 
quench h — )• —h. We will study the long-time behavior of both odd and even observables. 
In our numerical analysis, we have examined different values of the initial and final value of 
the parameter /i, and found that in the limits h <^ 1 and /i ^ 1, the quench dynamics are 
essentially trivial because the initial and final Hamiltonian share the same eigenvectors. For 
this reason, we will discuss only the intermediate case 

^prc— quench ^ y ^post— quench -j^ (4 3) 

Given these constraints we will still however consider two cases: 

• In one case we will look at ensembles of sparse random matrices, motivated by the 
previous considerations about the relationship between locality and sparseness. It should 
be stressed however that while a local observable will necessarily be sparse the converse is 
not necessarily true. Nevertheless the study of sparse random matrices may provide some 
reliable insights into some of the questions of the thermalization in local Hamiltonians. 

• In the second case we will look at matrices which are dense. Here then we are deviating 
from the motivation provided by the quantum Ising model. 

In both cases the entries of the random matrices will be generated according to the Gaussian 
orthogonal ensemble (GOE). 

We first consider quenches involving dense matrices. 



5 Thermalization in Dense Random Matrices Ensemble 

To define the Hamiltonian in the dense case, we generate three N/2 x N/2 matrices A,B,C 
and then assemble them according to Eq. (|4.1 p . The matrices A, C are symmetric and chosen 
according to the measure, fJ-{M), of a properly normalized GOE ensemble: 

f iVTr(M2)\ 
/i(Af)^exp(^ , 

while the matrix B has all of its entries distributed according to a normal distribution with 
mean and variance equal to For /i = ±1 the Hamiltonian itself will also be distributed 

according to the GOE ensemble and therefore the eigenvalues obey the semicircle law: 

p{E) = ^V^L^-E^ . (5.1) 

The spectrum thus falls in the range [— 2L,2L] and is therefore extensive as required. The 
observables are obtained with an analogous procedure, generating new matrices A,B,C and 
then using the expressions Eq. (|4.2 p . The numerical results reported below are calculated 
according to the following procedure: five instances of the Hamiltonian are generated according 
to the prescriptions above and for each instance of the Hamiltonian forty instances of the 
observable are generated. The relevant quantities are calculated for each instance of the 
observables, then the results are averaged. 
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Figure 5.1: Dense random matrices (with A'^ = 4000). Overlaps for the quench process. Left: 
I Cap for the 'ground' state. Right: |cap for the 2000*^^ state, in the middle of the energy band. 



5.1 Numerical results 

One of the prerequisites for the ETH to operate is that given the initial state I-^q) with energy 

_ 1 , , 

e = -^(V'ol-f^postlV'o) , 

the structure of its overlaps \ca\^ = KV'ol-^a)!^ with the post-quench eigenstates, as a function 
of the intensive energy — E^l is peaked around 6. 

We find this to be not the case for dense matrices, as can be explicitly seen by the two 
sample states in Fig. 15.11 drawn from the bottom and middle of the spectrum. Moreover, 
calculating the standard deviation of the energy on the initial state, we find that it is always 
large (around 1/4 of the range of the total spectrum) for all initial states, showing that the 
relation Eq. (|3.4 p does not hold, i.e. the effective dimension a satisfies a > 2. The broad 
distribution of overlaps is confirmed by the analysis of the Inverse Participation Ratio (IPR), 
defined as 

IPR = . (5.2) 

We show in Fig. 15.21 the IPR for the eigenstates of a single realization of a dense matrix. The 
IPR in this case is sharply distributed around N/2>. This finding can be understood through 
a simple model of random vectors on a N-sphere of unit radius (Porter- Thomas distribution). 
By a simple integration one finds |30| 



(5.3) 
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Figure 5.2: Dense random matrices (with N = 4000). Typical IPR of the initial states. 
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Figure 5.3: Dense random matrices. IPR vs. matrix size. The red lines are linear fits y = ax. 
Left panel: average IPR on all initial states (a = 0.3336). Right panel: maximum IPR 
(a = 0.3764). 



and therefore the IPR scales as 

^ N/3 . (5.4) 

This scaling is confirmed by our data, as shown in Fig. 15.31 Moreover, the fact that the mean 
IPR and the maximum IPR almost coincide is confirmation that all initial states are equivalent. 
This means that the pre-quench and the post-quench bases of the energy eigenvectors are 
completely random with respect one another. Eigenstates therefore have no reason to be 
localized in energy. 

Let us now turn our attention to the expectation values of observables since the main 
content of the ETH concerns the distribution of the eigenstate expectation values (EEVs), 
Oaa, and their behavior when the system size is increased. We first report two sample EEV 
distributions, given in Fig. 15. 4| which show no energy dependence. We can argue (and we have 
checked numerically) that the distribution of an intensive observable over the whole energy 
spectrum shrinks to zero for increasing system size. Moreover, it is not only the variance 
but even the support of the distribution of the observables that goes to zero inasmuch as the 
difference of the EEV maximum and minimum is going to zero as — )• oo (see Fig. 15.8 p . 

More precisely, we have: 

Oaa = {Ea\0\Ea) = P ' (5-5) 




-2-1012 -2-1012 

Figure 5.4: Dense random matrices. EEV {Ea\0\Ea) vs Ea- Left: even observable. Right: 
odd observable. 
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Figure 5.5: Dense random matrices. EEV variance (averaged over the entire spectrum) vs. 
matrix size. The continuous Hne is the fit a/x^ with b = 0.(9). The data points for the even 
and odd observables exactly overlap in this plot. 



where f3 indexes the eigenstates, of the observable O, while is the corresponding 
eigenvalue, and Aa^i3 = \{Ea\P)\'^. To estimate the r.h.s. we argue for an equivalency of 
observables and hold that the IPR of an eigenvector \Ea) of the post-quench Hamiltonian 
relative to the basis of eigenvectors |/3) equals the IPR of the initial state {tpo) in the basis 
\Ea). So we can suppose that Oaa can be expanded in terms of a set of ^ states, each of 
which is given by 

A. ,^^^r^^ . (5.6) 



" IPR{H O) 

Then, if we assume Aafi and Aap' are independent and note that Oaa has zero mean, we obtain 



\0 |2 ~ _ 



?,\N ) 



34 

N ' 



where cr^ is the variance of the spectrum of the observable O. In Fig. 15. 5| the numerical 
results are plotted together with a power-law fit and, as expected, the exponent is indeed 
close to one. 

Now let us consider the full support of the distribution of the EEVs where we define 60 as 
the difference of the maximum and the minimum of the EEVs among all the energy eigenstates 




1000 



2000 3000 
Matrix Size 



4000 



1000 



2000 3000 
Matrix Size 



4000 



Figure 5.6: Dense random matrices. Max-min EEV vs. matrix size for the range of even 
(left) and odd (right) observables (again the behavior of the two is indistinguishable). The 

continuous line is the fit a\ 



InN 
N ■ 



12 



7. X 10"^ 
5.x 10"^ 

3.x 10-" 

2. X lOi 
1.5 X 10"" 

1.x 10"*' 




1000 



1500 2000 
N 



3000 



7. 
5. 



. x 10"^ 
. x 10-*" 

3. x 10"^ 

% X IQ"^ 
1.5 X 10-*' 

1.x 10"" 



1000 



1500 2000 
N 



3000 



Figure 5.7: Dense random matrices, a = ^ (Omic 

laying in the central part of the spectrum e ~ 0. The continuous line is the fit a/x'' 
even observable h = 1.(9); right: odd observable b = 2.(1) 



vs. matrix size for initial states 

Left: 



\Ea)- Since the distribution is symmetric about zero, we have: 

6o = 2msix{Oaa} ■ (5.7) 

a 

To estimate the scaling of this quantity, we again approximate all the overlaps ^a,/3 as in Eq. 
(15.6 p . Therefore we are led to estimate the maximum of the quantity 

3 ' 



N 

where the prime on the sum indicates that only 1/3 of the total /3's are being summed over. 
As before we suppose that the random variables Ojs are independently distributed according 
to the intensive semicircle law (O^ G [—2,2]): 

p{x) = Prob(C'« = x) = ^y/A-x"^ . (5.8) 

In this case we can use large deviation theory (see for instance |31) ) and it follows that 

Prob(Oaa > x) ~ e-^^(^) , (5.9) 
where I{x) is the rate function given by 

I(x) = sup[9x- X{e)] . (5.10) 

e 

The function A(x) is the cumulant generating function 

e^(^) = £p(x)e^^ = ^ , (5.11) 

with Ii the modified Bessel function of first kind. Now the distribution governing the proba- 
bility that 5o is less than a value M is given by 

Prob(fo < M) = n ProhiOaa < M) ^ (l - e-^^(^^)) ^ 1 - ATe'^^^^) . 

a=l 
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Figure 6.1: Sparse random matrices (with N=4000). Left: density of states, right: level 
spacing statistics, the continuous hne is the Wigner surmise for GOE matrices. 



We can find the scaling of the typical value of the maximum by requiring that the probability 
is large enough 

IniV 

Prob((5ci < M) ~ const. =^ I{M) ~ . 

Since I{x) > 0, we are interested in the behavior of I{x) for small x and from Eq. (|5.10 P 

and Eq. (|5.11 P we get 

/in N 

l{x)^- + 0{x') 5o • (5.12) 

In Fig. 15. 6| one sees the numerical agreement with our heuristic argument. 

Finally let's consider the behaviour of the difference between the diagonal and the mi- 
crocanonical ensembles with increasing matrix size. To this end, we analyzed the difference 
a = \J (Omicro — Cdiag)^ of the two ensembles for each initial state Due to the broad en- 

ergy distribution of the overlap of the initial states, their intensive energy e = (■0o|-f^(— ^)|V'o) 
lies in a region, [—0.5,0.5], smaller than the range of the post-quench energies —2 < Ca < 2. 
a show the same behavior independent of the particular initial state I-^q) being considered, as 
can be argued by the constant IPR combined with a structureless EEV distribution. As all 
initial states are equivalent, we focus our attention on those initial states belonging to a small 
energy window around e = 0: the result is shown in Fig. 15.71 As can be seen, the difference 
between micro-canonical and diagonal ensemble rapidly goes to zero as a function of N. 

In conclusion, quenches in dense random matrices are characterized by initial states with 
large IPRs, EEV distributions of the post-quench eigenbasis with no energy dependence and 
whose variance goes to zero exponentially with increasing system size. In this sense, their ther- 
malization is trivial, as the spread of the micro-diagonal ensemble is governed by a distribution 
whose support is increasingly localized near zero as system size grows. 



6 Thermalization in Sparse Random Matrices 

We now turn to the more interesting case of thermalization in sparse random matrices. As we 
have indicated such matrices describe the Hamiltonians of systems with local interactions. In 
order to define the ensemble of these matrices, we employed a symmetric mask matrix, Jv[. 
This matrix has I's on the diagonal and in each of its rows we allow it to have on average IniV 
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Figure 6.2: Sparse random matrices. Behavior of the overlaps |cq 
(left) and for one in the upper portion of the spectrum (right). 



for a state midspectrum 



off-diagonal entries equal to 1. All remaining entries of the mask matrix are equal to zero. See 
Appendix B for details of how A4 is defined. The upper triangular part of the (symmetric) 
Hamiltonian is then obtained as: 

jj(h) _ / with di drawn from M{0,lnN) if i = j; . ^, 

V )t<j — I Q.. X ^ji-]^ ij.. cjrawn from Af{0, 1) if i < j, ^ ' ' 

where A/'(/x,cr^) is a normal distribution with mean and cr^ variance. Then the coefficients 
in the off-diagonal blocks are multiplied times h to reproduce the structure in Eq. (14.1 p . The 
different choice for the variances of the diagonal di and off-diagonal elements Oij is motivated 
by the requirement that the spectrum be extensive. In fact we can compute the variance of 
the spectrum: 

cx|,-^^^-21niV. (6.2) 

In the approximation in which the eigenvalues Ea are considered independent and normally 
distributed, one can relate the minimum of the spectrum with the variance, obtaining for the 
groundstate energy the estimate 

Egs = inin{Ea} ~ -anV^AnN ~ -21niV. 
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Figure 6.3: Sparse random matrices (with N = 4000). The IPR for a specific realization. The 
ordering of the initial states in this plot is according to their energy relative to the post-quench 
Hamiltonian, (^/^ol-f^'^lV'o)- 
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Figure 6.4: Sparse random matrices. IPR vs. matrix size. The continuous line is the y = ax fit. 
Left, average IPR on the whole spectrum (a = 0.0738). Right, maximum IPR (a = 0.2096). 



The eigenvalues midspectrum will be at most a few standard deviations an from the average 
0. Thus a typical eigenstate will be such that 

-4s^ oclniV . (6.3) 

We see then that our choices satisfy the requirement that energy is an extensive quantity. 

We compare these estimates with numerics in Fig. 16. II where we plot the density of states 
and the level spacing distribution for one realization of a sparse matrix. The density of states 
is no longer a semicircle, looking rather more like a bell-shaped distribution. Moreover, notice 
that, unlike the GOE case where the intensive quantities have a finite distribution in the large 
limit, i.e.Eq. (15.8 p . in this sparse case the (intensive) standard deviation behaves as ^/^^ , 
while the support of the spectrum remains approximately [—2,2]. We also see from the right 
side of Fig. 16.11 that the level-spacing distribution obeys the Wigner form for a non-integrable 
model. 

To generate the observables we follow a similar procedure, i.e. we employ the same mask 
M. The idea behind this choice is that the matrix M is responsible for the local structure in 
the Hilbert space and so we keep it for all the physical observables (including the Hamiltonian 
itself). So we have for the upper triangular part of a symmetric observable 

„ _ ( di with di drawn from Af{0, 1/ In N) i = j; 

~ \ Oij X Mij with Oij drawn from J\f{0, l/(ln N)"^) i / j. 
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Figure 6.5: Sparse random matrices (with N = 4000). The IPR of the post-quench eigenstates 
w.r.t. the local basis for a specific realization. The ordering of the initial states in this plot is 
according to their energy eigenvalue. 
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Figure 6.6: Sparse random matrices. IPR w.r.t. local basis vs. matrix size. The continuous 
line is the y = ax fit. Left, average IPR on the whole spectrum (a = 0.0321). Right, maximum 
IPR (a = 0.1632). 



The even and odd parts are then obtained as before by splitting O into diagonal and off- 
diagonal blocks. 



6.1 Numerical results 

Unlike the dense case, in sparse random matrices the overlap distributions are peaked around 
the initial state post-quench intensive energy e, as can be seen in two examples shown in Fig. 

The IPR is no longer constant (as it was for the random dense matrices), but shows 
behaviour dependent on the energy of the initial state, as demonstrated in Fig. 16.31 However 
there is still scaling of the IPR with the matrix size, as can be seen by studying the behaviour 
of the maximum IPR vs the matrix size, plotted in Fig. 16.41 In this case the mean and the 
maximum IPR are rather different, due to the presence of states with very small IPR. 

We see a similar phenomena when we study the IPR of the pre-quench and post-quench 
(they are statistically equivalent) eigenbasis relative to the local basis (that is, the basis of 
states in which the Hamiltonian matrices, H(h'^) and H{h^) are expressed). We see in Fig J6.6l 
that the mean and max IPR's here as a function of matrix size, A^, are similar to those in Fig. 
El 

The equivalency of these two different IPRs is not surprising. If we pick an eigenstate of the 
pre-quench Hamiltonian, |V'o)) that has a small IPR relative to the local basis, it will necessarily 



1.0 

A 0-5 
= 0.0 

0.5 

1.0 



V 



-2 



«• • • 





* Ti y^jt • • 

.'••.♦«•' 
*. " ■ • 


. ■ 







0.10 
0.05 
0.00 
-0.05 
-0.10 









• ■» :'* ' 







Figure 6.7: Sparse random matrices. EEV {Ea\0\Ea) vs Ca- Left: even observable. Right: 
odd observable 
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Figure 6.8: Sparse random matrices. EEV variance vs. matrix size for the full spectrum. 
The continuous lines are the fit ajx^. Left: even observable h = 0.5(0) Right: odd observable 
6 = 0.7(5). 



be only weakly coupled to the off-diagonal blocks, particularly as the Hamiltonian matrices are 
sparse, {i/jq) will then be closely related to some post-quench eigenstate. This implies in turn 
that {iPq) will have a small IPR relative to the post-quench eigenbasis. Similarly a pre-quench 
eigenstate |?/^o) with a large IPR in terms of the local basis will be strongly affected by the 
quench in the sense that it is unrelated to any eigenstate of the post-quench eigenbasis, and 
so will have a large IPR in terms of this basis. It is this that lies behind the similar shapes of 
Fig. 13] and Fig. ESI 

The behavior of the IPR relative to the local basis (and, by this equivalency, the IPR rel- 
ative to the post-quench eigenbasis) can be understood directly in the framework of Anderson 
localization |37| . Even though our problem is a many body one, our sparse Hamiltonian can be 
seen as akin to the dynamics of a non-interacting particle hopping on a Bethe-lattice of fixed 
connectivity, InA^, where each site has a random potential (the diagonal part of the random 
Hamiltonian). It is well known |36] that for this model the Anderson localization transition 
occurs with the presence of a mobility edge which separates the delocalized states (in the 
middle of the band) from the localized states (in the tails of the energy spectrum). In our 
case localization occurs when the IPR is 0(1), while delocalization is seen for IPR = 0{N). 

The position, Em, of the mobility edge can be determined by the following equation, 
derived along the same lines as |36] 

/■oo 

21niV / {p'{x-Em)-p'{Em-x)}lnxdx = l , (6.4) 
Jo 
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Figure 6.9: Sparse random matrices. Max-min EEV vs matrix size for the full spectrum. Left: 
even observable. Right: odd observable. 
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Figure 6.10: Sparse random matrices. EE V variance vs. matrix size for a small energy window 
around e = 0. The continuous line is the fit aj^ . Left: even observable h = 1.2(9). Right: 
odd observable b = 1.(2). 



where p{x) is the probability density of the diagonal terms in our Hamiltonian 

States with an \E\ > \Em\ such that the right hand side of Eq. (|6.4 p is less than 1 are 
localized. Otherwise they are delocalized. The integral in Eq. (16.4 p can be estimated at the 
leading order in the large limit, and one obtains 

(F^ \ / 9 1 n M 
InlniVW ^ E„ ~ iVlniVlnlniV . (6.6) 
2 m A* / V vr 

Thus in the large limit all the states with a non-zero intensive energy E/hiN behave as 
localized. Nonetheless the majority of the states, as they are concentrated in a window of 
width an (given in Eq. 16.2 p and as Em/cTH ^ 1, will be delocalized with an IPR = 0{N). 

The localization of states with finite intensive energies has implications for the EEV dis- 
tribution vs. Ca- As the observables have a matrix structure closely related to that of the 
Hamiltonians, we expect that localized states to be close to eigenstates of the observables 
itself, in contrast with Eq. (|5.6 P in the dense case. On such states, the observables will 
have expectation values far from their zero average. On the contrary, the delocalized states 
midspectrum will have EEV values closer to the mean of zero. In Fig. 16.71 we see numerical 
verification of this. 
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Figure 6.11: Sparse random matrices. Max-min EEV vs. matrix size for a small energy 
window around e = in logscale. The continuous line is the fit a/x^ + c. Left: even observable 
e = 0, & = 0.6(5), c = 0.06(7). Right: odd observable e = 0, 6 = 0.5(0), c = 0.0(1). 
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This result marks a strong difference witli respect to tlie case of dense matrices. We also 
see marked differences between the sparse and dense cases with both the variance and the 
support of the EEVs distribution of the observable, as shown in Fig. 16.81 and 16.91 while 
the variance approaches zero as N grows, the support does not, instead tending to a non-zero 
constant. Therefore the scaling Eq. (|5.12 p is no longer applicable, most likely as the overlaps, 
Aj^^, and the eigenvalues of the observable, O^, in Eq. (|5.5 P can no longer be considered as 
independent. 

Since the distribution of the overlap coefficients, c^, are peaked around the energy e = 
('0o|-f^(^)|'0o)j it is worthwhile to analyze the scaling behavior of the distribution of the EEVs 
in the vicinity of a specific e. Here we choose two different energy windows, one centered 
around e = 0, lying exactly mid-spectrum and one around e = 1. For e = the variance and 
the max-min spread as functions of the size N are plotted in Fig. (|6.10p and Fig. (|6.1ip . It 
is this distribution that is going to determine whether with respect to a particular observable 
we see thermalization. We again see that the variance is going to zero, while in contrast to 
the full spectrum, the max-min spread seems to tend, asymptotically, to either a very small 
constant value or to zero. The errors in our numerics are then not small enough to tell us at 
e = whether there is a complete absence of rare states. However we can be more definitive 
for the energy window centered at e = 1. In this window we see (as evidenced in Fig. I6.12P 
that the max-min spread of the EEVs in the large N limit tends to a finite constant for both 
the even and odd observables. 

Our numerical data then suggests that there are rare states where the observable remains 
far from its average value in each microcanonical energy window corresponding to non-zero 
intensive energy, differently for what has been observed in the t — J model |33| [3l| 135) . In 
these studies of the t-J model, rare states were argued to be absent for a range of strengths of 
the next nearest neighbour interaction, V', and for an energy window centered mid-spectrum. 
In particular with V small the t-J model is effectively integrable and rare states were found 
to be present and with V' large, the model develops energy bands, also compatible with the 
existence of rare states. It is for a middle range of V' that rare states were then found to be 
absent. Our own study sees some agreement with these results. And it is natural to think 
there would be some agreement at least. Our study of random matrix model should correspond 
to the t — J model with intermediate values of V': our matrix models are neither integrable 
nor do they have any notion of energy bands. In our case, rare states (at least for what we 
call the odd observable) seem to be vanishing as system size increases exactly in the center 
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Figure 6.12: Sparse random matrices. Max-min EEV vs. matrix size for a small energy 
window around e = 1. The continuous line is the fit a/x^ + c. Left: even observable e = 1, 
b = 0.3(5), c = 1.(7). Right: odd observable e = 1, 6 = 0.3(3), c = 0.1(4). 
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Figure 6.13: Sparse random matrices, a = (Omicro ~ Odiag)^ vs. matrix size for initial 
states laying in the central part of the spectrum e ~ 0. The continuous lines are the fits a/x^. 
Left: even observable h = 1.1(5). Right: odd observable h = 1.3(0). 
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Figure 6.14: Sparse random matrices, a 

states laying in a small window around e = 1. The continuous lines are the fits a/x^. Left: 
even observable b = 0.3(1). Right: odd observable b = 0.4(4). 



of the spectrum. However away from this midpoint of the spectrum, we find that rare states 
do exist, even in the thermodynamic limit. It would thus be interesting to extend the work of 
|33| IMl ES] to additional energy windows. 

Following |32], we then conclude in the case of sparse matrices that thermalization may 
depend on the particular nature of the initial state and will not occur when such rare states 
are given a proportionally large weight in the decomposition of the initial state. We, however, 
do not find for the particular initial conditions specified by our quench protocol that the rare 
states are given disproportional weight such that thermalization does not occur. For both 
energy windows e = (see Fig. (16.13P ) and e = 1 (Fig. (16.14P ). we see that with increasing 
matrix size the difference between the diagonal and microcanonical ensembles averaged over 
all initial conditions tends to zero. This implies that the weighting of rare states in our initial 
states is not preponderant. We do note however that the vanishing of the difference between 
ensembles decreases considerably more slowly with system size for the energy window, e = 1, 
than for e = 0. We might ascribe this to the presence of rare states at this energy - even 
though these states do not lead to non-thermalization in the thermodynamic limit, they may 
slow the approach to a thermalized state as the system size is increased. 

We verify this by computing the Kullback-Leibler entropy. This entropy is an information 
theoretic tool used to estimate how close two distributions are. It is defined as 

5kl = E ^(") 1^ , (6.7) 



21 



where P{a) is the expected distribution and Q{a) is the distribution to be compared. Skl is 
zero if the two distributions coincide except for sets of zero measures. In our case we choose 
P{a) = I Cap and Q{a) to belong to either a uniform or a Gaussian distribution centered about 
the energy e. The range over which Q{a) is defined has been taken such that the variances of 
Q and P coincides. Fig. 16.151 shows the average KL-entropy vs. matrix size for the central 
part of the spectrum which indicates in both cases a slow decay as the system size is increased. 
The slow approach of the distribution P{a) to a distribution Q{a) that is both smooth and 
symmetric about e suggests that rare states are not weighted in a peculiar way, thus permitting 
thermalization. 

In conclusion, for sparse random matrices our numerical data is compatible with the exis- 
tence of rare states. However, the initial states selected by the quench protocol do not seem 
to have large overlaps with these rare states and so we typically find thermalization as the 
end result of our quench process. 

7 Time Scales of Thermalization 

The sparse random ensemble, inasmuch as it mimics some characteristic properties of ther- 
malization in systems with local Hamiltonians, is the right framework to address the study 
of the thermalization time. Interest in this quantity can be traced back to the seminal paper 
by Von Neumann |6] regarding the quantum ergodic theorem (QET). The statement made in 
von Neumann's paper is that, under suitable assumptions (the Hamiltonian has no resonances 
- meaning that the energy level differences are non degenerate), any state \iPq) in the energy 
shell [e — A, e -|- A], will thermalize for most choices of the observable and most times t, i.e. 

Ot = (V'o(i)IClV'o(i)) = Tr(Opnic) for almost all t,0 . 

To make the notion of (macroscopic) observable and "most" used here precise would require 
the development of an involved technical apparatus and so instead we refer the reader to the 
existing literature [6l[7]. Nonetheless we can say there are important differences between this 
quantum thermalization and the classical notion of ergodicity where a time-average is involved. 
It is therefore of interest to give a precise estimate of the time needed for thermalization. There 
have been different approaches which have tried to clarify this question |39 | I40 | E | I42). In Fig. 
17.11 we examine the typical behavior of a realization of a random observable, Ot- It decays 
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Figure 6.15: Sparse random matrices. Kullback-Leibler entropy vs matrix size for the center 
of the spectrum for the uniform (squares) and the Gaussian distribution (circles). 
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Figure 7.1: Left: time evolution of the observable Oj versus time t. Right: time evolution of 
the fluctuations AP/Aq* versus time t, logscale on y axis. The straight hne is the exponential 
fit whose slope defines rg. 



towards the average value given by the diagonal ensemble: 

Ooo = Tr(Opdiag) , 

and we define the time Tp as the first time at which Ot meets Ooo- We stress here, that even 
if this quantity has not a direct physical interpretation, it can be considered as a lower bound 
for the thermalization time. The time evolution of the observable however keeps fluctuating 
around the average, due to finite system size, with it coming close to its initial value after a 
time. Tree, the recurrence time. To quantify these fluctuations we define 

Af ^ J [\o, - O^fdr ^ -y^Ol, = AS , 

a,o 

where we have assumed the absence of energy degeneracies and resonances. The quantity Af 
can be considered as the variance in the time interval [0, t\ of the observable and its behavior 
is plotted in Fig. 17.11 relaxation to the infinite time value is found as expected. We can fit 
this curve supposing an exponential relaxation, e"*/"^^ , defining in this way another time scale, 
Ts, the time interval needed for the relaxation of the fluctuations. Notice that this quantity is 
the one closer to the Von Neumann formulation: indeed, from the Chebyshev inequality, one 
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Figure 7.2: Scaling with the matrix size of the two time scale Tg (left) with a linear fit 
(-3.94 + 0.066A^) and Tp (right) with a logarithmic fit (-28.61 + 8.60 In A^). 
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has a bound on the fraction, /i/t, of times where the observable has an expectation vakie far 
from its average 

/i(re[0,t] and jO^ - Cool > a) Af 
t ^ ^ ■ 

As in the thermodynamic Umit — )• 0, this fraction must also go to zero in the long time 
limit. In Fig. 17. 2| we see the comparison between the two timescales ts and Tp versus the 
system size. From this plot one can see clearly that ts is a long time-scale, with a behavior 
proportional to N, i.e. the size of the Hilbert space. Notice that this can also be interpreted 
as the minimum spectral gap and at the leading order in A^: 

minl^;^ - Eb\ oc ^. 

In contrast, Tp is characterized by a much slower scaling with the size of the system and is 
therefore a fast time-scale. Although it is not easy to extract the precise scaling law from the 
available data, we have fit this data with the form 

Tp{N) = alnN , 

and so taking the scaling of this time scale to go as the volume. In contrast, the time scale 
for dense matrices was recently argued to go as the inverse of the volumejH]. 



8 Conclusions 

In this paper we have addressed the issues of thermalization and the Eigenstate Thermaliza- 
tion Hypothesis in the framework of random matrices, aiming to identify certain statistical 
properties of quantum extended systems subjected to a quench process. For this purpose we 
focused our attention on Z2 breaking quantum Hamiltonians, among the simplest theoretical 
quench protocols. In an attempt to encode in our analysis the property of locality, we have 
considered the ensemble of sparse random matrices and we have compared the data coming 
from this ensemble with similar data extracted from the ensemble of dense random matrices. 
We have found reliable evidence of different behavior in the two ensembles. These differences 
show up both in the IPR of the quench states and in the distribution of the expectation values 
of the observables on post-quench energy eigenstates. In particular, while in the dense random 
matrix ensemble both the variance and the support of the observables vanish with increasing 
system size L, the sparse random matrix ensemble sees instead strong indications that the 
variance of EEVs goes to zero while the support remains finite as L — )■ 00. The different 
behavior of the two ensembles can be traced back to the different density of states exhibited 
by the two sets of matrices: while in the dense matrices all states are delocalized in the Hilbert 
space, with almost equal overlap on all energy eigenstates, in the sparse matrices there are 
instead both delocalized and localized states. Localized states give rise to rare values of the 
expectation values of the observables, i.e. values which differ from the typical ones sampled 
by the micro-canonical ensemble. If properly weighted, such localized states may give rise to 
a breaking of thermalization. In the absence of such weighting, as seems to be the case in 
the initial conditions chosen by our quench protocol, one instead observes relaxation to the 
thermal value of the local observables. 

In the framework of sparse random matrices, we have also provided numerical estimates 
of the different time scales of thermalization. We have found that it is possible to identify two 
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time scales: a fast one Tp and a slow one ts, and that they depend differently upon the size 
of the system. 
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A Total momentum and quasiparticles 



We want to consider in this appendix the effect of a change of basis on the sparseness of the 
Hamiltonian matrix. In particular we will focus on a change to a basis where momentum is 
a good quantum number. We take once again the Ising guiding example. Changing 

the basis, the new matrix representation of the Hamiltonian Eq. (|2.1 p can be obtained in 
terms of an unitary matrix U : 

n' = u^nu (1.1) 

and, for a general change of basis, the resulting matrix will not be sparse anymore; this is quite 
intuitive: for an arbitrary unitary matrix U, i.e. an arbitrary basis, any notion of locality is 
lost. 

However we can ask what happens in a much more common case: is the Hamiltonian 
still sparse in the basis of momenta? We will show that there is a subtle issue related to the 
definition of this change of basis. To be more specific, let's put the Ising Hamiltonian on a 
circle by introducing periodic boundary conditions such that the (L + l)-th site is identified 
with the first one. In this case the system becomes translationally invariant; this can be 
formally written as: 

[H,T]=0, 

where T is the shift operator, defined by: 

T\mi, . . .,mL) = \mL,mi, . . .,mL-i). 

Clearly we have: 

= 1, (1.2) 
and therefore we can define the total momentum P as: 

T = e'^ ^ [H, P] = 0. 

As T is a unitary operator, using Eq. (|1.2 p we easily deduce the usual structure of the 
spectrum of P as The basis of momenta is then the basis of the eigenstates of the total 
momentum. However, the eigenvalues of P are highly degenerate and so many definitions 
are possible. Two examples will demonstrate that sparseness depends on which definition is 
employed. 

Let us first construct a complete set of eigenstates of P in the following way. We first pick 
a state \s) = \mi, . . . mi) in the real-space basis. We obtain an invariant subspace for T by 
considering the set of states: 

= Span{|s),r|s),r2|s),...,r^»-i|s)} 
where Lg is been defined asH 

Lg = minjre such that T'^\s) = |s)} 
The eigenstates of P composed of linear combination of states in the set Is can be written as: 

e—T'ls) P|s„)= |s„). (1.3) 



"^The minimum exists because the set contains at least the element, L. We also note that Ls divides L. 
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A full basis of eigenstates for P can be obtained repeating this procedure for different states 
\s): we call this basis the rigid-translation Fourier basis (RTFB). It is easy to understand that 
the resulting matrix in this new basis is still sparse. In fact each state is a superposition of at 
most L states and it follows the corresponding transformation U contains at most L non-zero 
entries in each row and column. Using Eq. (jl.l p . we conclude that the Hamiltonian matrix 
in this basis is still sparse having at most non-zero entries in each row (L^ ^2''"). 

Is this the end of the story? Let us consider what happens if we use a different momentum 
basis. To this end we work in a second-quantization framework and focus upon a set of 
operators satisfying (where the it stands for the fermionic and bosonic case): 

[ai,aj]± = Sij. 

These operators create and destroy a single particle excitation at position i. A real-space basis 
can be written in this formalism as: 

{a\r{alr-..\n). (1.4) 

The same formalism can be adopted in the Ising case by setting^ at ~ and taking 
= I 4 . . . 4). Here we can define an excitation with defined momentum by setting: 



L-l 
j=0 



Putting the inverse of this expression into Eq. (|1.4 p . we get a new basis of eigenstates of total 
momentum P: 

J J io\ - J J 




We call this the single-particle Fourier basis (SPFB). However once we restrict to a subspace 
where P is defined (appearing as a block for the Hamiltonian matrix), there is a strong 
difference between this case and the one defined in Eq. (|1.4 p : in fact here, not only the full 
state but each component excitation has a defined momentum. To better understand this 
difference we consider a two particle case (in the one particle case, there are no differences 
between the two momentum bases). Consider a state \s) = SxiSx2\^) = \xi,X2), that is the 
state with only two up spins in positions, xi and X2- In the RTFB, this state can be expressed 
as a linear combination of L states, \sn), with n = 0, . . . , L — 1, where |s„) are obtained as 
superpositions of the rigid translations, T^\s). In these translated states the distance \x2 — xi\ 
between up states always remains the same. On the other hand, with the SPFB, each two- 
particle state has a non-zero matrix element with \xi,X2)- 

Increasing the number of particles to M, in the RTFB case we always have one summation 
with ~ L terms. Instead in the SPFB, we will have M summations, corresponding to ~ L^^ 
terms. 

The conclusions we can draw from these considerations are that: 



*In Irf this transformation can be made more rigorous using the Jordan-Wigner transformation. 
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• A local Hamiltonian will appear as a sparse matrix in the real-space basis. 

• If we consider the Fourier basis obtained as a superposition of the rigid translations of the 
real-space basis, the Hamiltonian will appear again as a sparse matrix though perhaps 
with a larger density of non-zero entries. This is true simply because the change of basis 
we are considering is sparse. 

• If we consider the Fourier basis obtained taking the Fourier transform of single particle 
states, the change of basis in each block of defined total momentum will not be sparse. 
In the general case, the Hamiltonian matrix will be characterized by dense blocks of 
fixed total momentum. 



B Generation of the mask 

Our aim is to generate a mask, telling us where the non-zero entries are, according to the 
scaling Eq. (|2.4 p . The most obvious choice would be simply to take the mask such that 
each non-diagonal entry is chosen independently, being one with probability p = and zero 
otherwise: Pq = 1 — p. But a simple computation shows that with this choice the probability 
Pnc of having a column (or row) with all the non-diagonal entries equal to zero will be quite 
high 

To avoid this, we try to fix the number of non-zero entries in each row. However, since at 
the end, we will need a symmetric matrix with one on the diagonal, we follow a slightly 
more complicated procedure. We first generate "half" of the mask matrix and at the end we 
symmetrize by summing it with its transpose. To be more precise, for each row, i, of the 
matrix we generate a set, composed of integers drawn from the set {1, . . . , z, . . . , A^}, i.e. a 
set where the element i is missing. The size, n^, of Ri has probability Pj„ of equaling n where 
we define by 

Pin = {q]K~l,[q\ + (1 - {<l]¥nM\ 

where {q} and [(^J represent respectively the fractional and integer part of q and q is defined 
by: 

q = [N -I) - yJ{N - 1)2 - {N -l)\nN. 

In this way, the size of the Ris will be either [(/J or \_q\ -|- 1, but such that the average is 
fii = q. The sets i?j's contain the indexes of the non-zero elements in each row. We now define 
the mask joining each row with the corresponding column, thus symmetrizing the matrix: 




1 {i = j) V {j e R-i) V {i e Rj); 
otherwise. 



Our choice of q takes into account the possibility of elements overlapping between the row and 
the corresponding column, ensuring we have the correct number of non-zero elements in each 
row: ^ 

#{i I Mij + 0} = 1 + 2(/ - -^—^ = 1 + IniV. 

Notice that the number of non-zero elements in each row is not fixed, because small fluctuations 
are induced by the symmetrization procedure; however we can be sure that it will be greater 
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than [q\ . We recall that In N/N is exactly the threshold for the connectedness of the Erdos- 
Renyi graph |38) : with this procedure there will be no isolated subblocks, as is confirmed by 
the Wigner-Dyson level statistics shown in Fig. 16.11 
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